##Intragroup_Conflict_Leadership_Study_3.R
##Study 3 - opposition research and impeachment (GOP Sample)

##Load data
setwd("C:/Users/lmh735/Dropbox (Political Science)/Partisanship and leadership-Trump/2021 POBE Submission/R&R/Replication Files")
data3<- read.csv("Study_3.csv", header=TRUE)

##set directory for output
setwd("C:/Users/lmh735/Dropbox (Political Science)/Partisanship and leadership-Trump/2021 POBE Submission/R&R/Replication Files/Output")

#########################################
##Appendix Table F2
##Balance using chi-squared, only Republicans
library(MASS) ##(only run for this, causes issues with tdyr/dplyr, need to detach if)
##gender
balance.gender<- table(data3$female, data3$treatment.rebuke.gop)
prop.table(balance.gender, 2)
chisq.test(balance.gender)
##education
balance.educ<- table(data3$educ, data3$treatment.rebuke.gop)
prop.table(balance.educ, 2)
chisq.test(balance.educ)
##income
balance.income<- table(data3$income, data3$treatment.rebuke.gop)
prop.table(balance.income, 2)
chisq.test(balance.income)
##age (categorical)
balance.age.cat<- table(data3$age.cat, data3$treatment.rebuke.gop)
prop.table(balance.age.cat, 2)
chisq.test(balance.age.cat)

##sample demographics
##gender
sample.gender<- table(data3$female)
prop.table(sample.gender)
##education
sample.educ<- table(data3$educ)
prop.table(sample.educ)
##income
sample.income<- table(data3$income)
prop.table(sample.income)
##age (categorical)
sample.age.cat<- table(data3$age.cat)
prop.table(sample.age.cat)

detach("package:MASS", unload=TRUE)

######################################
##Appendicx Table F3
#summary stats
summary(data3$trump_ft_01)
length(na.omit(data3$trump_ft_01))
sd(data3$trump_ft_01, na.rm=TRUE)

summary(data3$senator_ft_01)
length(na.omit(data3$senator_ft_01))
sd(data3$senator_ft_01, na.rm=TRUE)

summary(data3$dissent_appr)
length(na.omit(data3$dissent_appr))
sd(data3$dissent_appr, na.rm=TRUE)

summary(data3$voteconroy)
length(na.omit(data3$voteconroy))
sd(data3$voteconroy, na.rm=TRUE)

summary(data3$loyalty)
length(na.omit(data3$loyalty))
sd(data3$loyalty, na.rm=TRUE)


#######################################
##Libraries for analysis
#Analysis with dplyr
library(dplyr)
library(tidyr)
library(infer)
library(ggplot2)
library(gridExtra)
library(tidyverse)

#######################################
##Appendix Table F1
##N by condition
impeach.N<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, trump_ft_01) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(n=n())
impeach.N
sum(impeach.N$n)

#####################################
##trump FT, post-Senator response
impeach_trump_ft_01<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, trump_ft_01) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(mean=mean(value),
            sd=sd(value),
            n=n(),
            se=sd/sqrt(n))

impeach_trump_ft_01

##plot
fig7a<- ggplot(impeach_trump_ft_01) + # define aes for each individual component instead
  geom_errorbar(aes(x=ResponseCondition.new, ymin=mean-1.96*se, ymax=mean+1.96*se), ##95% CI
                width=.15) +
  geom_point(aes(x=ResponseCondition.new, y=mean), size=2) + # add the points on
  #facet_grid(~ Condition.short) + # split into plots, side by side
  scale_y_continuous(limits=c(0,1)) + 
  theme(axis.text=element_text(size=12), axis.title.y = element_text(size = 12),
        plot.title = element_text(size=12))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ylab("Feeling Thermometer") +
  ggtitle("(a) Trump Feeling Thermometer")

##t-test
t.test(as.numeric(data3$trump_ft_01) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")

######################
##Senator Conroy FT
impeach_senator_ft_01<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, senator_ft_01) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(mean=mean(value),
            sd=sd(value),
            n=n(),
            se=sd/sqrt(n))

impeach_senator_ft_01

##plot
fig7b<- ggplot(impeach_senator_ft_01) + # define aes for each individual component instead
  geom_errorbar(aes(x=ResponseCondition.new, ymin=mean-1.96*se, ymax=mean+1.96*se), ##95% CI
                width=.15) +
  geom_point(aes(x=ResponseCondition.new, y=mean), size=2) + # add the points on
  #facet_grid(~ Condition.short) + # split into plots, side by side
  scale_y_continuous(limits=c(0,1)) + 
  theme(axis.text=element_text(size=12), axis.title.y = element_text(size = 12),
        plot.title = element_text(size=12))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ylab("Feeling Thermometer") +
  ggtitle("(b) Senator Feeling Thermometer")

##t-test
t.test(as.numeric(data3$senator_ft_01) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")

##########################
##Approval of Senator's position
impeach_dissent_appr<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, dissent_appr) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(mean=mean(value),
            sd=sd(value),
            n=n(),
            se=sd/sqrt(n))

impeach_dissent_appr

##plot
fig7c<- ggplot(impeach_dissent_appr) + # define aes for each individual component instead
  geom_errorbar(aes(x=ResponseCondition.new, ymin=mean-1.96*se, ymax=mean+1.96*se), ##95% CI
                width=.15) +
  geom_point(aes(x=ResponseCondition.new, y=mean), size=2) + # add the points on
  #facet_grid(~ Condition.short) + # split into plots, side by side
  scale_y_continuous(limits=c(0,1)) + 
  theme(axis.text=element_text(size=12), axis.title.y = element_text(size = 12),
        plot.title = element_text(size=12))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ylab("Approval") +
  ggtitle("(c) Approval of Senator's Statement")

##t-test
t.test(as.numeric(data3$dissent_appr) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")

#######################
##Primary vote for Senator
impeach_voteconroy<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, voteconroy) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(mean=mean(value),
            sd=sd(value),
            n=n(),
            se=sd/sqrt(n))

impeach_voteconroy

##plot
fig7d<- ggplot(impeach_voteconroy) + # define aes for each individual component instead
  geom_errorbar(aes(x=ResponseCondition.new, ymin=mean-1.96*se, ymax=mean+1.96*se), ##95% CI
                width=.15) +
  geom_point(aes(x=ResponseCondition.new, y=mean), size=2) + # add the points on
  #facet_grid(~ Condition.short) + # split into plots, side by side
  scale_y_continuous(limits=c(0,1)) + 
  theme(axis.text=element_text(size=12), axis.title.y = element_text(size = 12),
        plot.title = element_text(size=12))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ylab("Likelihood of Vote for Senator") +
  ggtitle("(d) Vote for Senator in Primary")

##t-test
t.test(as.numeric(data3$voteconroy) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")

#############################
##loyalty index (mechanism)
impeach_loyalty<- data3 %>%
  mutate(ResponseCondition.new=factor(treatment.rebuke.gop, levels=c("Loyalist", "Dissent"), labels=c("Loyalist", "Dissent"), ordered=FALSE)) %>%
  select(ResponseCondition.new, loyalty) %>%
  gather("measure", "value", -ResponseCondition.new) %>%
  filter(!is.na(value)) %>%
  group_by(ResponseCondition.new) %>%
  summarize(mean=mean(value),
            sd=sd(value),
            n=n(),
            se=sd/sqrt(n))

impeach_loyalty

##plot
fig7e<- ggplot(impeach_loyalty) + # define aes for each individual component instead
  geom_errorbar(aes(x=ResponseCondition.new, ymin=mean-1.96*se, ymax=mean+1.96*se), ##95% CI
                width=.15) +
  geom_point(aes(x=ResponseCondition.new, y=mean), size=2) + # add the points on
  #facet_grid(~ Condition.short) + # split into plots, side by side
  scale_y_continuous(limits=c(0,1)) + 
  theme(axis.text=element_text(size=12), axis.title.y = element_text(size = 12),
        plot.title = element_text(size=12))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ylab("Loyalty Index") +
  ggtitle("(e) Senator Loyalty Index")

##t-test
t.test(as.numeric(data3$loyalty) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")

##########################################
##output all DVs together
grid.arrange(fig7a, fig7b, fig7c, fig7d, fig7e, nrow=2)
fig7.g<- arrangeGrob(fig7a, fig7b,fig7c, fig7d, fig7e, nrow=2)
ggsave("Fig7_impeach_all_DV.png", fig7.g, height=10, width=15)
##Impeachment (Jan 2020)

##########################################
##Appendix Table I1
##Power calcs for Senator response

##Senator FT
##loyal
summary(data3$senator_ft_01[data3$treatment.rebuke.gop=="Loyalist"])
sd(data3$senator_ft_01[data3$treatment.rebuke.gop=="Loyalist"])
length(data3$senator_ft_01[data3$treatment.rebuke.gop=="Loyalist"])
##Dissent
summary(data3$senator_ft_01[data3$treatment.rebuke.gop=="Dissent"])
sd(data3$senator_ft_01[data3$treatment.rebuke.gop=="Dissent"])
length(data3$senator_ft_01[data3$treatment.rebuke.gop=="Dissent"])
##Power analysis
library(stats)
power.t.test(n=89, sig.level=0.05, power=0.80, sd=0.28)
##Note used lower of N in loyal and N in dissent for n, sd is larger sd in loyal and dissent
##could have detected a mean difference of 0.12 (observed difference was 0.27)
t.test(as.numeric(data3$senator_ft_01) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")
0.3077528 - 0.5791089 

##Senator vote
##loyal
summary(data3$voteconroy[data3$treatment.rebuke.gop=="Loyalist"])
sd(data3$voteconroy[data3$treatment.rebuke.gop=="Loyalist"])
length(data3$voteconroy[data3$treatment.rebuke.gop=="Loyalist"])
##Dissent
summary(data3$voteconroy[data3$treatment.rebuke.gop=="Dissent"])
sd(data3$voteconroy[data3$treatment.rebuke.gop=="Dissent"])
length(data3$voteconroy[data3$treatment.rebuke.gop=="Dissent"])
##Power analysis
library(stats)
power.t.test(n=89, sig.level=0.05, power=0.80, sd=0.34)
##Note used lower of N in loyal and N in dissent for n, sd is larger sd in loyal and dissent
##could have detected a mean difference of 0.14 (observed difference was 0.307)
t.test(as.numeric(data3$voteconroy) ~ data3$treatment.rebuke.gop, na.rm=T, alternative="two.sided")
0.3632697 - 0.6699356 

##########################################
##Footnote
##Comparison of pre-treatment FT to ANES
summary(data3$trump_ft_pre)
sd(data3$trump_ft_pre)
